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SUMMARY 

The NASA STRu ctural Analysis (NASTRAN) (Ref 1) program is one of the most extensively 
used engineering applications software in the world. It contains a wealth of matrix operations and 
numerical solution techniques, and they were used to construct efficient eigenvalue routines. The 
purpose of this paper is to examine the current eigenvalue routines in NASTRAN and to make 
efficiency comparisons with a more recent implementation of the Block Lanczos algorithm by 
Boeing Computer Services (BCS). This eigenvalue routine is now available in the BCS mathe- 
matics library as well as in several commercial versions of NASTRAN. In addition, CRAY main- 
tains a modified version of this routine on their network. Several example problems, with a 
varying number of degrees of freedom, were selected primarily for efficiency bench-marking. 
Accuracy is not an issue, because they all gave comparable results. The Block Lanczos algorithm 
was found to be extremely efficient, in particular, for very large size problems. 

INTRODUCTION 

In NASTRAN the real eigenvalue analysis module is used to obtain structural vibration modes 
from the symmetric mass and stiffness matrices, and K/^, which are generated in the pro- 
gram using finite element models. Currently the user has a choice of four methods for solving 
vibration mode problems: Determinant Method, Inverse Power Method with Shifts, Tridiagonal 
Method (Givens’ Method) and Tridiagonal Reduction or FEER Method. NASTRAN provides all 
these options for user convenience as well as for analysis efficiency. For example, the Givens’ 
Method is most appropriate when all the eigenvalues are of equal interest. By the same token, it is 
not suitable (because of the need for excessive computer resources) when the number of degrees 
of freedom is too large (greater than three to four hundred) unless preceded by Guyan reduction 
(ASET or OMIT). The Inverse Power, Determinant and FEER Methods are most suitable when 
only a small subset of the eigenvalues are of interest. These methods take advantage of the 
sparseness of the mass and stiffness matrices and extract one or a small subset of eigenvalues at a 
time. 
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•NASTRAN without qualification refers to COSMIC-NASTRAN (or government version) in the paper. 
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The purpose of this paper is to examine, in some detail, the real eigenvalue analysis methods cur- 
rently available in NASTRAN and to make efficiency comparisons with the Block Lanczos algo- 
rithm as implemented by Boeing Computer Services (BCS) and currently available in some 
commercial versions of NASTRAN (for example MSC-NASTRAN and UAI-NASTRAN). The 
accuracy of the eigenvalues is not an issue in this paper, because all the methods gave comparable 
results Efficiency in terms of computer time is the only issue in this bench-marking. This study 
was made, for all cases, on a single platform, the CRAY XMP. The genesis of the Block Lanczos 
Method in all the NAS TRANs, as well as the CRAY version, is the one implemented by BCS with 

some modifications. 

Section 1 discusses the general form of the eigenvalue problem for vibration modes. In Section 2 
a mathematical formulation of the four methods in NASTRAN is given with emphasis on the 
FEER Method as a precursor to the Lanczos Method. A detailed mathematical description of the 
Block Lanczos Method is given in Section 3. Also reference is made to the Lanczos method in 
MSC NASTRAN and to its implementation by CRAY Research, Inc. In Section 4 selected fre- 
quencies are calculated for five structures of varying complexity using the Inverse Power Method, 
the FEER Method, MSC/NASTRAN Lanczos Method and CRAY Lanczos Method. Results are 
discussed in Section 5 and recommendations are made for possible implementation into NAS- 
TRAN. 
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1.0 The Eigenvalue Problem 

1.1 The general form of the eigenvalue problem for vibration modes is 


Kx = XMx 


( 1 ) 


2 

where M and AT are the symmetric mass and stiffness matrices, the eigenvalue X = © the square 
of the natural vibration frequency, and * is the eigenvector corresponding to X. The dimension of 
the matrices AT and Af is nxn, where n is the number of degrees of freedom in the analysis set. For 
this paper it is assumed that AT and Af are at least positive semi-definite. Thus associated with Eq 
(1) are n eigenpairs X., such that 

i (2) 


(3) 


where Af,-, is referred to as the modal mass or generalized mass. It is evident from Eq (3) that the 
eigenvectors are orthonormal with respect to the mass matrix. Also the eigenvectors are orthonor- 
mal with respect to the stiffness matrix, i.e. 


Kx } = 

X.Mx i i- 1,2 

Properties of the eigenvectors include: 
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xJKxj = 


K H ioT i= J 


(4) 


v 0 for i * j 

where AT,, is the modal stiffness or generalized stiffness. 

The Rayleigh quotient shows that the modal mass, Af„, and modal stiffness, AT,-,-, are related to the 
eigenvalue X., i.e. 


X. = 

l 


x T i Kx i 

x^MXj 


K ii 

M : 


(5) 


ll 


For normalized eigenvectors with respect to modal mass, Eqs (3) can be written as 

f 1 for i =J 


x]Mx j = 


Now using Eqs (5), Eqs (4) can be written as 


xJKxj = 


0 for /' * j 


X . for = j 
J 

0 for / * j 


( 6 ) 


(7) 
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The central issue of a real eigenvalue or normal modes analysis is to determine the eigenvalues, 
and the eigenvectors, x i7 which satisfy the conditions stated by Eqs (1-7). The next sections 
present the important elements of the eigenvalue methods of interest. 
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2.0 Eigenvalue Extraction Methods in NASTRAN 

2.1 For real symmetric matrices there are four methods of eigenvalue extraction available in 

NASTRAN: the Determinant Method, the Inverse Power Method with shifts, the Givens’ 

Method of Tridiagonalization and the Tridiagonal Reduction or FEER Method. Most methods of 
algebraic eigenvalue extraction can be categorized as belonging to one or the other of two groups: 
transformation methods and tracking methods. In a transformation method the two matrices M 
and K are simultaneously subjected to a series of transformations with the object of reducing them 
to a special form (diagonal or triadiagonal) from which eigenvalues can be easily extracted. 
These transformations involve pre and post multiplication by elementary matrices to annihilate 
the off-diagonal elements in the two matrices. This process preserves die original eigenvalues in 
tact in the transformed matrices. The ratio of the diagonal elements in the two matrices gives the 
eigenvalues. In a tracking method the roots are extracted, one at a time, by iterative procedures 
applied to the dynamic matrix consisting of the original mass and stiffness matrices. In 
NASTRAN the Givens’ and the FEER methods are transformation methods, while the 
Determinant and the Inverse Power methods are tracking methods. Both tracking methods and 
the Givens’ method will be discussed briefly in this section while the Lanczos algorithm, the main 
emphasis of this paper, is outlined here and in more detail in the next section. 

2.2 Determinant Method 


For the vibration problem 


Kx = XMx 


( 8 ) 


the matrix of coefficients, A, has the form 

A = K-XM (9) 

The determinant of A can be expressed as a function of X, i.e. 

D(A) = \A\ = (X-X { ) (X-X 2 ) ...(X-XJ 

where X., i = 1,2 ...n are the eigenvalues of A. In the determinant method D(A) is evaluated for 
trial values of X, selected according to an iterative procedure, and a criterion is established to 
determine when D(A) is sufficiently small or when X is sufficiently close to an eigenvalue. The 
procedure used for evaluating D(A) employs the triangular decomposition 

A = LU ( 10 ) 
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for an assumed value of X where Lisa lower unit triangular matrix and U is an upper triangular 
matrix. D(A) is equal to the product of the diagonal terms of U. Once an approximate eigenvalue, 
X,., has been accepted, an eigenvector, is determined from 

LUx i = 0 ( u > 

by back substitution where one of the elements of x, is preset. Since I (X,) is nonsingular, only 
U (X.) is needed. The determinant method may not be efficient in some cases if more than a few 
eigenvalues are desired because of the large number of triangular decompositions of A. 

2.3 Inverse Power Method with Shifts 

The Inverse Power Method with shifts is an iterative procedure applied direcdy to Eq (1) in the 
form 

[K-XM]x = 0 (12) 


It is required to find all the eigenvalues and eigenvectors within a specified range of X. Let 


X = X +A 
o 


(13) 


where \ is a constant called the shift point. Therefore A replaces X as the eigenvalue. The iter- 
ation algorithm is defined in the nth iteration step by: 


[K-\M]w n = Mx n _ 1 


1 


X y. = — W 

n r ft 
n 


(14) 

(15) 


where c„, a scaler, is equal to that element of the vector w n with the largest absolute value. At 
convergence l/c„ converges to A, the shifted eigenvalue closest to the shift point, and x„ con- 
verges to the corresponding eigenvector <j> r Note from Eq (14) that a triangular decomposition of 
matrix K - XAf is necessary in order to evaluate w„. The shift point X o can be changed in order to 
improve the rate of convergence toward a particular eigenvalue or to improve accuracy and con- 
vergence rates after several roots have been extracted from a given shift point. Also \ a can be 
calculated such that the eigenvalues within a desired frequency band can be found and not just 
those that have the smallest absolute value. 

For calculating additional eigenvalues, the trial vectors. .v„, in Eq (14) must be swept to eliminate 
contributions due to previously found eigenvalues that are closer to the shift point than the current 
eigenvalue. An algorithm to accomplish this is given as follows: 


m 


*n = 5,- I <*>*»>*< 


(16) 


i = 1 
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where x n is the trial vector being swept, m is the number of previously extracted eigenvalues, and 
<{>. is defined by 



X U N 

x J,N Mx i,N 


(17) 


where Xj N is the last eigenvector found in iterating for the ith eigenvalue. 

The inverse power method allows the user to define a range of interest on the total fre- 

quency spectrum and to request a desired number of eigenvalues, ND, within that range. When 
ND is greater than the actual number of eigenvalues in the range, then the method guarantees the 
lowest eigenvalues in the range. 

2.4 Givens’ Method of Tridiagonalization 


In the Givens’ method the vibration problem as posed by Eq (8) is first transformed to the form 


Ax = Xx 


(18) 


by the following procedures. The mass matrix, M, is decomposed into upper and lower triangular 
matrices such that 



(19) 


If M is not positive definite, the decomposition in Eq (19) is not possible. For example, when a 
lumped mass model is used, NASTRAN does not compute rotary inertia effects. This means that 
the rows and columns of the mass matrix corresponding to the rotational degrees of freedom are 
zero resulting in a singular mass matrix. In this case the mass matrix must be modified to elimi- 
nate the massless degrees of freedom. 

Thus Eq (8) becomes 

Kx = XLL T x (20) 


/ T 1 

which implies after premulitplying by L and post multiplying by ( L ) that 

L ~ 1 K (L T ) \ = Xx 


( 21 ) 


i.e. 


Ax = Xx 


where A-L' 1 K(L* Y 1 . Note that L 1 is easy to perform, since L is triangular. Also 

— I T ~ ^ 

A = L K{L ) is a symmetric matrix. The matrix A is then transformed to a tridiagonal 


148 


matrix, A n by the Givens’ method, i.e a sequence of orthogonal transformations, 7y, are defined 
such that 

T r T r -l-‘ T 2 T l Ax = kT r T r _ 1 ...T 2 T l x (22) 

Recall that an orthogonal transformation is one whose matrix T satisfies 

TT 7 = 7 T T = / ( 23 ) 


the identity matrix. The eigenvalues of A arc preserved by the transformation, and if 


x 


rpT 'T'T tT 

1 l i 2*" i r- 



(24) 


then from Eq (22) 


T r T r -l- T 2 T l AT I T 2- T *-l T ry = XT r T r - V - T 2 T JVt ' ^ 


i.e. 

T r T r _ v ..T 2 T l AT\T T 2 ...T r r _ ] T T r y = Xy (25) 

by repeatedly applying Eq (23). Eq (25) implies that y is an eigenvector of the transformed matrix 
T r T r _ v ..T 2 T i AT r l Tl...T r r _ l T T r . Thus x can be obtained from y by Eq (24). 

The eigenvalues of the tridiagonal matrix, A p are extracted using a modified Q-R algorithm, i.e., 
A r+ j = Q T r A r Q r such that A r is factored into the product QfR r where R r is an upper triangular 
matrix and Q r is orthogonal. Thus 

A r = Q r R r (26) 


A r+ 1 _ 2 r ^ r G r 

- Q T r Q r R r Q r 

Since Q r is orthogonal, then 

A r + 1 = 


from Eq (26) 


(27) 


In the limit a sr->» and A is symmetric, A r will approach a diagonal matrix. Since eigenvalues 
are preserved under an orthogonal transformation, the diagonal elements of the limiting diagonal 
matrix will be the eigenvalues of the original matrix A. 

To obtain the ith eigenvector, yj, of the tridiagonal matrix, A p the tridiagonal matrix A r — \J is 
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factored such that 


(28) 


A r - U = L-Uj 

where I, is a unit triangular mania and U, is an upper triangular mania. The eigenvector y, is then 
obtained by iterating on 

U.yW »»> 

where the elements of the vector v are arbitrary. Note that the solution of Eq (29) is easily 
obtained by back substitution since £/, has the form 

Pill r l 

p 2 q 2 r 2 

v,= - - - 

Pn-l^n- 
p n 

The eigenvectors of the original coefficient matrix, A, are then obtained from Eq (24). 

Note that in the Givens’ method the dimension of A equals the dimension of A r The major share 
of the total effort expended in this method is in converting A to A r Therefore the total effort is not 
strongly dependent on the number of eigenvalues extracted. 

2.5 Tridiagonal Reduction or FEER Method 

The tridiagonal Reduction or FEER method is a matrix reduction scheme whereby the eigenval- 
ues in the neighborhood of a specified point, \ o , in the eigenspectrum can be accurately deter- 
mined from a tridiagonal eigenvalue problem whose dimension or order is much lower than that 
of the full problem. The order of the reduced problem, m, is never greater than 

m = 2q + 10 

where q is the desired number of eigenvalues. So the power of the FEER method lies in the fact 
that the size of the reduced problem is the same order of magnitude as the number of desired 
roots, even though the actual finite element model may have thousands of degrees of freedom. 

There are five basic step in the FEER method: 

1. Eq (8) is converted to a symmetric inverse form 

Bx = A Mx < 31 > 



150 


where 


A = ^ (32) 

o 

and A, is a shift value. 

O 

2. The tridiagonal reduction algorithm or Lanczos algorithm is used to transform Eq (31) into a 
tridiagonal form of reduced order. 

3. The eigenvalues of the reduced matrix are extracted using a Q-R algorithm similar to that 
described in Section 2.4. 

4. Upper and lower bounds on the extracted eigenvalues are obtained. 

5. The corresponding eigenvectors are computed and converted to physical form. 

To implement Step 1, consider Eq (8), 

Kx = XMx 

When vibration modes are requested in the neighborhood of a specified frequency, X Q , Eq (8) can 
be written 

Kx — X Mx — XMx — X Mx 

o o 

(K - X q M) x = ( X-X q )Mx (33) 

Let K = K- X q M and X' - X - X Q . Then from Eq (33) 

Kx = X'Mx (34) 

_-l 

x = X'K Mx 

_-l 

Mx = X'MK Mx 

MK l Mx = Mx (35) 

Factor K by Cholesky decomposition, i.e. 

K = Ld'lJ (36) 

where L is a lower triangular matrix and d' is a diagonal matrix. Then Eq (35) can be written 


-1 



Bx = AMx 


where B = m[ (L t ) V 1 !" 1 ]*/ and A = ^ Now step 1 is complete. 

O 

To implement Step 2 rewrite Eq (31) as 

Bx = Ax 


where B = M~ l B . Now B is reduced to tridiagonal form, A, using single vector Lanczos recur- 
rence formulas defined by 


a U = 


yjBVf 


Vi+i = BV r«u v r d i v i - 1 

1/2 


\i = 1, 2, m 


(37) 


—T 

■ t r* 


{v; + iMv f+ i> 


V 


1+ 1 




i = 1, 2, ...m - 1 


where vector V o =0, V; is a random starting vector and dj=0. The reduced tridiagonal eigenvalue 
problem is now given as 


Ay 


a ll d 2 

d 2 a 22 

d 3 a 33 d 4 

\ \ \ 

d m - 1 a m — 1, m — 1 d m 

d a 

m mm 


Ay 


(38) 


where A approximates the eigenvalue A of Eq (31), and y is an eigenvector of A. The Lanczos 
formulas generate a V matrix, vector by vector, i.e. 

V = [V lf V 2f ...VJ (39) 

and Eqs (37) are modified by NASTRAN such that each vector V i+1 is re-orthogonalized to all 
previously computed V vectors, i.e. V is orthonormal to M. 

V T MV = / (40) 
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Thus 


A = V T BV 


(41) 


Note from Eq (41) that A is an mxm matrix. 

For step 3 the eigenvalues. A, and eigenvectors, y, of Eq (38) are obtained as described for the 
Givens’ method in Section 2.4. The eigenvectors are normalized so that 

yfy,- = 1 1 = 1 , ...,m (42) 


For step 4 the following error bound formula has been derived and serves as a criterion for select- 
ing acceptable eigensolutions 


> 4 - 

< 

+1 ymi 

a. 

i 


A,(1+VV 


(43) 


In Eq (43) X. is an approximation to the exact eigenvalue X ( . in Eq (8), d m + x is calculated from 
Eqs (37), y mi is the last component of the mth eigenvector, y m , of A, and A t is the ith eigenvalue 
of A. The ith eigenvalue X ( is acceptable, if £. is less than or equal to a preset error tolerance. 


Now step 5 is implemented for acceptable eigenvalues. If (A, y) is an eigenpair of Eq (38), then 

Ay = Ay 


or from Eqs (40) and (41) 

V T BVy = A V T MVy 


BVy = A MVy 


(44) 


Now if x =VV, then 


Bx = AMx 


i.e. (A, x) is an eigenpair of Eq (3 1 ). 

Thus for step 5 the eigenvectors of Eq (3 1 ) or equivalently Eq (8) are calculated from 

-v = Vy (45) 


and the eigenvalue X is calculated from Eq (32) i.e. 

X = I+X 
A ° 


(46) 


Note that in the FEER method the matrix B enters the recurrence formulas, Eqs (37), only through 
the matrix -vector multiply terms BV r Therefore B is not modified by the computations. Lanczos 
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procedures for real symmetric matrices require only that a user provide a subroutine which for 
any given vector, z, computes Bz. 
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3.0 Block Lanczos Method 


3.1 Recall that the eigenvalue problem in vibration analysis is given by Eq (8), i.e. 

K.x = XMx 

where K and M are symmetric positive definite matrices. Generally the eigenvalues of interest are 
the smallest ones, but they are often poorly separated. However, the largest eigenvalues which are 
not interesting have good separation. Also convergence rates are very slow at the low end of the 
spectrum and fast at the higher end. Convergence rates can be accelerated to the desired set of 
eigenvalues by a spectral transformation, i.e. consider the problem 

M(K — gM)~ 1 M.x - uMx (47) 

where a, the shift, is a real parameter. It can be shown that ( X, x) is an eigenpair of Eq (8) if and 
only if (r-^— ,.r) is an eigenpair of Eq (47). The spectral transformation does not change the 
eigenvectors, but the eigenvalues of Eq (47) are related to the eigenvalues of Eq (8) by 

" ■ xh (48) 

This transformation will allow the Lanczos algorithm to be applied even when M is semidefinite. 

Consider the effect of the spectral transformation on a satellite problem which will be discussed in 
detail in Section 4. Figure 1 shows the shape of the transformation. Table A shows the effect of 
the transformation using an initial shift of c = .046037. Note that the smallest 8 eigenvalues are 
transformed from closely spaced eigenvalues to eigenvalues with good separation. 
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Problem 
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.07229 

38.09088 

.03611 
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22.05574 
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1 
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4 ; 

.31296 

3.74640 
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.00084 

2.291 14 x 10' 3 
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.40134 

1.88521 

.05142 

6 

.58357 

| 1.86035 

j .16180 

.24002 

i .43042 

j .01174 

“T 

/ 

.74537 

: 1 .42993 

.00103 

.00153 

.00210 

5.72784 x 10° 


S .74640 1 .42783 

Table A 

Our objective is to define the Spectral Transformation Block Lanczos algorithm. Let s consider 
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first the Basic Block Lanezos Algorithm. 

3.2 Basic Block Lanezos Algorithm 

Consider the Lanezos Algorithm (Refs 2.3) for the eigenvalue problem. 

Hx = \x (49) 


where H is symmetric 

The block Lanezos iteration with block size p for an nxn matrix H is given as: 

Initialization: 
set Q 0 = 0 
set Bj = 0 

choose Rj and orthonormalize the columns of Rj to obtain Qj 
Lanezos Loop: 

For j — 1,2,3, ... 
set Uj = HQj- Qj.j Bj 
set Aj = Qj Uj 

set R j+1 = u r G/Aj 

Compute the orthogonal factorization Qj+jBj + j = Rj+j 
End Loop 

Matrices Qj, Uj. and Rj for j- 1, 2 , ... are nxp; Aj and Bj arc pxp. Aj is symmetric and Bj is upper 
triangular. The blocksize p is the number of column vectors of Qj. So if p = 1 , then Qj is a column 
vector, q. Thus the matrix H is not explicidy required, but only a subroutine that computes Hq for 
a given vector q. Aj and Bj are generalizations of the scalers aj and dj in the ordinary Lanezos 
recurrence. 

The recurrence formula in the Lanezos loop can also be written as 

R j+ 1 = Qj+ 1 R j + 1 = HQj- QjAj-Qj_ ( 50 ) 

The orthogonal factorization of the residual. implies that the columns of Qj are orthonormal. 
Indeed it has been shown that the combined column vectors of the matrices. Qj. Q 2 . ... Qj, called 
the Lanezos vectors, form an orthonormal set. 

The blocks of Lanezos vectors form an nxjp matrix Wj where 

Wj = [Q v Q 2 Qj) ( 51 ) 
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From the algorithm itself a jpxjp block tridiagonal matrix, Tj, is defined such that 



^1*2 

0 ... o 

^2^2 

#3 ... 0 

0 ... 

1 

i 

0 ... 

0 Bj Aj\ 


(52) 


Since the matrices B } are upper triangular, Tj is a band matrix with half band width p+1 . The first 
j formulas defined by Eq (50) can be combined using Eqs (51) and (52) into a single formula 

HWj - WjTj + Qj +1 Bj +1 eJ (53) 

where Ej is an nxp matrix of zeros except the last pxp block is a pxp identity matrix. 
Premulitplying Eq (53) by wj implies 

wjHWj = wJw j T J + wjQ j+1 B j+l Ej 

i.e. 

WjHWj = Tj (54) 

since 

wJWj = I and W jQj+i = 0 

Eq (54) implies that Tj is the orthogonal projection of H onto the subspace spanned by the col- 
umns of Wj. Also if (0, s) is an eigenpair of T y i.e. TjS=sB , then (X, W jS) is an approximate 
eigenpair of H. A discussion on the accuracy of the approximation will be delayed until the spec- 
tral transformation Block Lanczos Algorithm is considered. Basically the Lanczos algorithm 
replaces a large and difficult eigenvalue problem involving H by a small and easy eigenvalue 
problem involving the block tridiagonal matrix Tj. 

3.3 Spectral Transformation Block Lanczos Algorithm 

Since our primary consideration is vibration problems, consider the eigenprobem posed by 
Eq (47) i.e. 

M{K- GM) ~ 1 Mx = uMx 


158 



The Lanczos recurrence with block size p for solving Eq (47) is given by 
Initializ ation 
set Q 0 = 0 
set Bj = 0 

choose Rj and orthonormalize the columns of Rj to obtain Qj with Q^MQ j = I p 

Lanczos Loop 

For; = 1,2,3, ... 

set Uj = (K-oM)' 1 (MQ) - Qj.j Bj 
set Aj = U] ( MQj ) 
set R j+l = Uj-QjAj 
Compute Q j+l and (MQ J+l ) such that 

a) (2y+l^>+l = ^y +1 

b) = I p 

End Loop 

Note that the algorithm as written requires only one multiplication by M per step and no factoriza- 
tion of M is required. The matrices Qj are now M orthogonal, rather than orthogonal, i.e. 

QjMQj = / ( 55 ) 

Also the Lanczos vectors arc M orthogonal, i.e. 

wJmWj = i 

The recurrence formula in the Lanczos loop can also be written as 

Q J+t B j+l = (K-cM)-'MQ ] -Q J A r Q j _ l Bj (56) 

Now, as before, combining all j formulas of Eq (56) into one equation yields 

(K-oM)-'MWj = WjTj + Qj +l B j+1 Ej (57) 

where Wj, Tj, and Ej arc as defined in Eq (53). Premulitplying Eq (57) by wJm implies 

wJm(K-<5M)~'mw j = wJmWjT + wjQj +l B j+l E] 

i.e. 

= Tj (58) 


159 


since 


wjMWj = / and Wj<2 j+l = 0 

Eq (58) implies that T, is M-orihogomd projection of (AT-oM)' 1 onto th. subspace spanned 

by the columns of IP,. The eigenvalues of 7} will approximate the eigenvalues of Eq (47). If 
(«, s) is an eigenpair of 7), then (0, Wfi will be an approximate eigenpatr of Eq (47). 

Recall that our main interest is in solving Eq (8). From Eq (48) 


or V = G + ^ (59) 

i.e. if 0 is an approximate eigenvalue of T Jt then from Eq (59) v is an approximate eigenvalue of 
Eq (8). Recall that the spectral transformation does not change the eigenvectors, therefore 
y = W s is an approximate eigenvector for Eq (8). 

Let's examine the approximations obtained by solving the block tridiagonal eigenvalue problem 
involving the matrix T } . Let (0, s) be an eigenpair of T j i.e. 

Tjs = sQ 

and let y = WjS. Then Premulitplying Eq (57) by M and post multiplying by s gives 
M (K - oM)~ i MWjS — MWjTjS = MQ j+1 B j+l Ejs 
M(K-aM)~ X My-MWjsB = MQ j+l B J+l Ejs 

M(K- - MyQ = MQ j+ ^B j+l E T j s (60) 

7 1 

Recall for any vector q, ll^ll^i = q ^ (Ref 4). 

Therefore, taking the norm of Eq (60) and using Eq (55) 

\\M(K-cM)~ l My-MyQ\\^ t = \\MQj +l Bj +l EjS ||^, 

= HVi E J 5 ^j (61) 

Note that p is easily computed for each eigenvector 5. It is just the norm of the p vector obtained 
by multiplying the upper triangular matrix B ■ + j with the last p components of s. 

From Ref 5 the error in eigenvalue approximations for the generalized eigenproblem is given by 


160 


|| M(K-oM) l My-MyQ || , 

1 < Ml = ft. < 62 > 

X=o " WMyW^ 

Thus is a bound on how well an eigenvalue of Tj approximates an eigenvalue of Eq (47). 

Recall that if 0 is an approximate eigenvalue of Tj, then from Eq (48) 

1 

V " G+ e 

is an approximate eigenvalue of Eq (8). Consider 

|X-V| = 

= s|< X - o) <jlo- e) 

(«) 

Therefore \\ - v| & . Thus 1 is a bound on how well the eigenvalues of Eq (47) approximate the 

e 2 e 2 

eigenvalues of Eq (8). 

3.4 An Analysis of the Block Tridiagonal Matrix Tj 

The eigenproblem for T- is solved by reducing Tj to a tridiagonal form and then applying the 
tridiagonal Ql algorithm. The eigenextraction is accomplished in three steps. 

1 An orthogonal matrix Q T is found so that T } is reduced to a tridiagonal matrix H , i.e. 

Q t t TjQt = H «*> 

2. An orthogonal matrix Q H is found so that H is reduced to a diagonal matrix of eigenvalues, 
A, i.e. 


Qh h Qh = A 


3. Combining Eqs (64) and (65) gives 

{QyfQj) T j (QjQtf) = A (66) 

where Q-j-Qh ls the eigenvector matrix for T } . The orthogonal matrices Q H and Q T are a product 
of simplex orthogonal matrices, Givens’ rotations, Qh^Qh^-Qh, or ^8°* 
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rithms used for steps (1) and (2) are standard and numerically stable algorithms drawn from the 
EISPACK collection of eigenvalue routines. 

Note from Eq (61) that only the bottom p entries of the eigenvectors of Tj are needed for the eval- 
uation of the residual bound. Therefore it is unnecessary to compute and store the whole eigen- 
vector matrix for T } . Only the last p components of the eigenvector matrix are computed. 

The error bounds on the eigenvalues Eq (62) and (63) are used to determine which eigenvectors 
are accurate enough to be computed. At the conclusion of the Lanczos run the EISPACK subrou- 
tines are used to obtain the full eigenvectors of Tj. Then the eigenvectors for Eq (47) are found 
through the transformation 

y = Wjs 

3.5 Other Considerations in Implementating the Lanczos Algorithm. 

The use of the block Lanczos algorithm in the context of the spectral transformation necessitates 
careful attention to a series of details: 

a. The implications of M-orthogonality of the blocks 

b. Block generalization of single vector orthogonalization schemes 

c. The effect of the spectral transformation on orthogonality loss 

d. The interactions between the Lanczos algorithm and the shifting strategy. 

All of these issues are addressed in detail in Refs. 5,6. 

3.6 The Block Lanczos algorithm as described in the previous sections was developed as a 
general purpose eigensolver for MSC NASTRAN (Ref 7). Boeing designed the software such 
that the eigensolver was independent of the form of the sparse matrix operations required to 
represent the matrices involved and their spectral transformations. The key operations needed 
were matrix-block products, triangular block solves and sparse factorizations. These, and the data 
structures representing the matrices, are isolated from the eigensolver. Therefore, the eigensolver 
code could be incorporated in different environments. 

For this paper we tested the block Lanczos algorithm as incorporated in MSC NASTRAN and as 
further developed by Boeing and incorporated into code by Cray Research, Inc. The block Lanc- 
zos algorithm in MSC uses the factorization and solve modules which are standard operations in 
MSC. The Cray Lanczos code uses the Boeing eigensolver with matrix factorization, triangular 
solves, and matrix -vector products from the mathematical libraries supplied by Boeing computer 
services (BSCLIB-EXT). For vibration problems the CRAY code can be used with the stiffness 
and mass matrices, K and M, as generated by NASTRAN. NASTRAN is run to generate binary 
files containing the K and M matrices. These files are input files to the Cray code which calculates 
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l ,* checks the orthogonality of the eigenvectors, at, via a: Kx, calculates the Rayleigh 
e.genvafoes^he ^ ^ eigenvalues. and calculates the no™ of the 

quotient at Kx/xMx P ^ ^ cigenvector files output from the CRAY 

eigenvector testdu . ^ ^ pr0 « ssing jf desired. Since the commerctal 

are suitable for tnpu NASTRANS do not give M and AT in the same formats, 

(MSC) and the government COSMIC) N R NASTRAN was used to repre- 

they need to be reformaned before calling the CRAY code. CSAR-N ASTRAN 

sent NASTRAN on the CRAY XMP. 
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4.0 Test Problems 


In this section several test problems were solved using the inverse power and FEER eigenvalue 
extraction methods in COSMIC NASTRAN, the Lanczos algorithm in MSC NASTRAN and the 
Lanczos algorithm as implemented by CRAY Research. These problems were chosen based on 
the complexity of the finite element model in terms of the kinds of elements used and the number 
of degrees of freedom. All methods as expected gave approximately the same numerical results. 
The only criterion used to compare the different methods was the number of seconds needed to 
reach a solution given that all problems were solved on the same platform, a CRAY XMP. 

4.1 Problem 1 Square Plate 

A square 200 in x 200 in plate in the x-y plane was modeled with QUAD4 elements only. Five 
meshes were defined. Details are given in Table 1. All elements were 1.0 in thick. Material prop- 
erties were constant for all meshes. Each plate was completely fixed along the x-axis and the y- 
axis at x=200 in. 


Number of Grid 
Points 

Number of 
Elements 

Number-of 
Degrees of 
Freedom 


MESH 


lOx 10 

20x20 

30x30 

40x40 

50x50 

121 

441 

961 

1681 

2601 

100 

400 

900 

1600 

2500 

515 

2015 

4515 

8015 

12515 


Table 1: DETAILS OF THE FIVE MESHES DEFINED ON THE SQUARE PLATE 

For all cases 5 frequencies were requested in the interval [0, 20hz]. Table 2 gives the results for 
the 10 x 10 plate, and Table 3 gives the results for the 50 x 50 plate. As expected within each case 
the numerical results from the different eigenextraction techniques are approximately the same. 
The differences in numerical results between the 10x10 case and the 50 x 50 case reflect the fine- 
ness of the mesh for the 50 x 50 case. Both Lanczos algorithms were ran with a fixed block size 
of p = 7. 
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r 



FREQUENCIES IN Hz 



1 

2 

3 

4 

5 

COSMIC Inverse 
Power 

6.2980 

7.1720 

11.6374 

17.4440 

18.3096 

COSMIC FEER 

6.2980 

7.1720 

11.6374 

17.4440 

18.3096 

MSC Lanczos 

6.2730 

7.2173 

11.7181 

17.2125 

18.3392 

CRAY Lanczos 

6.2730 

7.2173 

11.7181 

17.2125 

18.3392 


Table 2: 10 x 10 SQUARE PLATE 


FREQUENCIES IN Hz 



1 

2 

3 

4 

5 

COSMIC Inverse 
Power 

6.4048 

7.6103 

12.5487 

17.6764 

19.3642 

COSMIC FEER 

6.4048 

7.6103 

12.5487 

17.6764 

19.3642 

MSC Lanczos 

6.4054 

7.6159 

12.5599 

17.6745 

19.3739 

CRAY Lanczos 

6.4054 

7.6159 

12.5599 

17.6745 

19.3739 


Table 3: 50 x 50 SQUARE PLATE 

Table 4 gives the CPU time in seconds from the CRAY XMP needed to extract five frequencies 
for each case. Recall that the CRAY Lanczos algorithm needs to obtain the mass and stiffness 
matrices in binary form from NASTRAN. Thus the time given for this algorithm is the total time 
from two computer runs, i.e. the time to obtain the mass and stiffness matrices plus the time to ran 
the Lanczos algorithm separately. 


165 




CPU TIME IN SECS 


MESH SIZE 


COSMIC Inverse 
Power 

COSMIC FEER 
MSC Lanczos 
CRAY Lanczos 


10 X 10 

20x20 

30x30 

40 x 40 

50x50 

14.734 

50.936 

97.801 

197.769 

328.830 

8.085 

19.363 

39.877 

77.994 

132.179 

4.783 

13.641 

30.973 

59.283 

103.188 

4.174 

11.170 

23.785 

45.433 

78.009 


Table 4: CPU TIME IN SECONDS TO OBTAIN 5 FREQUENCIES 

Figure 2 is a plot of the degrees of freedom versus the CPU time in seconds on the CRAY for the 
four eigenvalue extraction techniques. 



Figure 2: Degrees of Freedom versus CPU Tune in Seconds. 
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4.2 Problem 2 Intermediate Complexity Wing 


A three spar wing shown in Figure 3 was modeled w«h 88 grids end 158 elements of *e Mow- 
ing typer 62 QUAD4, 55 SHEAR, 39 ROD and 2 TRIA3. All elements vaned m thickness 
J^crional area. Material pmperiies were the same for all elements. The *mg was com- 
pletely fitted at the root which left 390 degrees of freedom. Ftve frequencies were requested tn the 
interval [0, 300hz], Table 5 gives d« frequencies calculated and the CPU time tn seconds for the 
four eigenettttaction algorithms. As for Problem 1 bod, Lancaos algorithms were run wtth . fitted 

block size of p = 7. 



COSMIC 
Inverse Power 

COSMIC FEER 

MSC Lanczos 


CRAY Lanczos 


Figure 3: Intermediate Complexity Wing 
FREQUENCIES IN Hz 


CPU TIME IN 
SECONDS 


1 

2 

3 

4 

5 


46.574 

135.924 

176.813 

205.030 

254.713 

10.314 

46.574 

135.924 ! 

176.813 

205.030 

254.713 

8.085 

46.573 

135.918 

176.811 

205.029 

254.690 

4.886 

| 46.573 

135.918 

176.811 

205.029 

254.690 

4.873 


Table 5: INTERMEDIATE COMPLEXITY WING RESULTS 
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4.3 Problem 3 Radome 


A composite radome shown in Figure 4 was modeled with 346 grids and 630 elements of the fol- 
lowing types: 54 TRIA2, 284 BAR and 292 QUAD4. The QUAD4’s were both isotropic and 
composite with 46 elements isotropic and 246 elements modeled as four cross-ply unsymmetric 
laminates of 40, 38, 36, and 32 layers, respectively. Hie radome was completely fixed at the base 
which left 1782 active degrees of freedom. Ten frequencies were requested in the interval 
[0,100 hz]. Table 7 gives the frequencies calculated and the CPU time in seconds for the four 
eigenextraction algorithms. Both Lanczos algoirthms were run with a fixed blocksize of p = 7. 



FREQUENCIES IN Hz 


CPU 
TIME IN 
SECS 


COSMIC 
Inverse Power 
COSMIC 
FEER 
MSC 
Lanczos 
CRAY 
Lanczos 


1 

2 

3 

4 ! 5 

6 

7 

8 

9 

10 


56.325 

67.946 

69.290 

i 

81.486j90.835 

90.971 

92.074 

92.410 

93.365 

101.441 

63.986 

i 

'56.325 

67.946 

69.29018 1 .486(90.835 

90.971 

92.074 

92.410 

93.365 

101.441 

21.318 

1 

56.068 

66.958l68.2l3i80.843 89.7 15l90.248i90.768 

i 

91. 676|92.365 

i 

j 98.729 

17.768 

i 

56.068 

66.958 

68.213 

80.843189.715 

90.248 

(90.768 

91.676 

92.365 

98.729 

13.854 


Table 6: Radome Results 


168 



4.4 Problem 4 Satellite 


A satellite shown in Figure 5 was modeled with 2295 grids and 1900 elements distributed as 
shown in Table 7. 



Figure 5: Satellite 


ELEMENT TYPE 


Number of 
Elements 


ROD 

BEAM 

ELAS1 

ELAS2 

TRIA3 

QUAD4 

BAR 

HEXA 

PENTA 

RBE2 

15 

134 

30 

8 

45 

111 

297 

40 

56 

498 


Table 7: Satellite Element Distribution 


Sixteen different matenals were referenced, and 34 coordinate systems were used. All elements 
varied in thickness and cross-sectional area, and concentrated masses were added to selected 
grids. The satellite has 5422 active degrees of freedom. Fifty frequencies were requested in the 
interval [0, 20hz]. Table 8 gives every fifth frequency calculated and the CPU time in seconds for 
the four eigenextraction algorithms. Again both Lanczos algorithms were run with a fixed block 
size of p = 7. 
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FREQUENCIES IN Hz 


CPU 
TIME 
IN SEC 


COSMIC 

Inverse 

Power 

COSMIC 

FEER 

MSC 

Lanczos 

CRAY 

Lanczos 


1 

5 

10 

15 

20 

25 

30 

35 

40 

45 

50 






NO S< 

DLUTI 

ON IN 

2000 S] 

ECS 




.072 

1 

.313 

1.497 

1.663 

2.419 

5.414 

9.000 

10.974 

13.328 

17.474 

19.758 

294.759 

.072 

.313 

1.497 

1.634 

2.406 

5.417 

9.056 

10.975 

13.267 

17.104 

19.649 

121.065 

.072 

.313 

1.497 

1.635 

2.406 

5.418 

9.056 

10.975 

13.268 

17.111 

19.650 

81.016 


Table 8: SATELLITE RESULTS 
4.5 Problem 5 Forward Fuselage - FS 360.0 - 620.0 

A section of a Forward Fuselage from FS 360.0 to 620.0 shown in Figure 6 was modeled with 
1038 grids and 3047 elements distributed as shown in Table 9. 

Eleven different materials were referenced. All elements varied in thickness or cross-sectional 
area. The fuselage was fixed in the 123 directions at FS 620.0. The model had 6045 active 
degrees of freedom . Sixty frequencies were requested in the interval [0, 20hz]. Table 10 gives 
every fifth frequency calculated plus the last one and the CPU time in seconds for the four eigen- 
extraction algorithms. Both Lanczos algorithms were ran with a fixed block size of p = 7. 
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Figure 6: Forward Fuselage 


ELEMENT TYPE 

CONROD SHEAR TRIA3 QUAD4 
Number of Elements 1141 | 885 395 15 572 

Table 9: Forward Fuselage Element Distribution 


BEAM 
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FREQUENCIES IN Hz 


ll 



2 

















5.0 Summary and Recommendations 

The current real eigenvalue analysis capability in NASTRAN in quite extensive and adequate for 
small and medium size problems. In particular the FEER Method s performance is reasonable at 
least for the problems tested in this paper. However, the Block Lanczos Method as implemented 
by BCS is more efficient for all the problems. 

An analysis of Section 4 results clearly shows that the Block Lanczos Algorithm merits consider- 
ation for possible implementation into NASTRAN. Comparing CPU secs Table 4 implies that the 
CRAY Lanczos method runs 94% to 64% faster than the FEER method. Similarly from Tables 5, 
6, 8 and 10 the CRAY Lanczos runs 66%, 54%, 260% and 177%, respectively, faster than the 
FEER method. 

The comparisons are not near as striking when we consider the CRAY Lanczos and the MSC 
Lanczos. Comparing CPU seconds the CRAY Lanczos runs from .2% faster in Table 5 to 105.7% 
faster in Table 10. The difference in CPU time reported for these two methods can be attributed to 
two factors: (1) algorithm enhancements and (2) the Boeing Extended Mathematical Subprogram 
Library (BCSLIB-EXT) versus the standard mathematical modules in MSC. The CRAY Lanczos 
is based on [Ref 5] which is, most recent, dated July 1991. The MSC Lanczos is based on [Ref 6] 
which is dated 1986 plus subsequent updates by MSC. All problems were run under MSC NAS- 
TRAN Version 66a. Recent communications with Roger G. Grimes at Boeing, one of the devel- 
opers of the-shifted Block Lanczos algorithm, reveals that the Lanczos algorithm is continuously 
being refined and improved. 

The problems chosen to test the four eigenextraction methods while diverse in terms of the num- 
ber of degrees of freedom and element distribution were stable with no clusters of multiple eigen- 
values. The multiple eigenvalue problem and its relation to the user chosen blocksize, p, is 
discussed in detail in [Ref 5]. The authors conclude that based on timing results for the selected 
problems, the shifted Block Lanczos Algorithm should be considered for possible implementation 
into NASTRAN. 

Boeing Computer Services is reluctant to sell or lease their Block Lanczos routine to public 
domain programs such as COSMIC-NASTRAN or ASTROS. In view of this the authors recom- 
mend the following alternatives: 

• Modify the FEER Method from a single vector Lanczos algorithm to a Block 
Lanczos algorithm. 

• Obtain the Block Lanczos algorithm from an alternate source. 

• Provide links for calling subroutines from the commercial math libraries such as 
the BCS or CRAY library. 
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